/*! 
 ***************************************************************************
 * \file singlePhase2D.cpp
 * Single phase in porous media in 1D.
 *
 * \par Description
 * This code solves the 2D single phase equation for pressure:
 * \f$ \frac{\phi \mu c}{k}\frac{\partial p}{\partial t} = \nabla^2 p \f$.
 * This equation is solved in \f$ (x, y) \in [0,100cm] \times [0, 10cm] \f$.
 * The data used are: 
 * The initial condition for the temperature is \f$ P = 0 \f$
 * and the boundary conditions:  \f$ P_L=2 atm \f$ and \f$ P_R=1 atm\f$.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./singlePhase  \endverbatim
 ***************************************************************************
 *  \author    Luis M. de la Cruz [ Mon Mar 17 08:47:37 GMT 2008 ]
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,2>::huge ScalarField2D;

void exactSolution(ScalarField2D &T, double hx, double dx, double PL, double PR)
{
  int NI = T.ubound(firstDim);
  int NJ = T.ubound(secondDim);
  double coef = dx / hx;

  for(int i = 0; i <= NI; ++i)
    for(int j = 0; j <= NJ; ++j)
        T(i,j) = PL + (PR - PL) * i * coef;
}

int main() 
{
  double length_x, length_y, dt, tolerance = 1e-7; 
  int num_nodes_x, num_nodes_y, frequency, tdma_iter = 2000;
  double porosity, viscosity, permeability, compresibility;
  double PL, PR, Tmax;

// Reading data from the input file
    std::ifstream input_cin ("inputSingle2D"); 
    input_cin >> length_x >> length_y
	      >> num_nodes_x >> num_nodes_y
	      >> porosity
	      >> permeability
	      >> viscosity
	      >> compresibility
	      >> PL
	      >> PR
	      >> dt
	      >> Tmax
	      >> frequency;
    input_cin.close();

    Tmax += dt;

    double Gamma = permeability / (porosity * viscosity * compresibility);

    std::cout << "\n Gamma = " << Gamma;

    StructuredMesh<Uniform<double, 2> > mesh(length_x, num_nodes_x,
					     length_y, num_nodes_y);
    double dx = mesh.getDelta(X), dy = mesh.getDelta(Y);
    mesh.print();
    ScalarField2D p ( mesh.getExtentVolumes() );
    ScalarField2D pa( mesh.getExtentNodes()   ); 
    ScalarField2D pn( mesh.getExtentNodes()   );


    pn = 1.0; // initial condition on nodes
    pa = 1.0; // initial condition on nodes
    p = 1.0; // initial condition on volumes

    exactSolution(pn, length_x, dx, PL, PR);
    NumUtils::interpolateToVolumes(pa, pn);

    Range all = Range::all();
    p(0, all) = 2.0;
    pn(0, all) = 2.0; // boundary condition

    InOut::writeToFile_DX(pn, 0, "./DataSingle2D/pres.", dx, dy);
       
    SparseMatrix< Diagonal< double, 2> > A(num_nodes_x, num_nodes_y); 
    ScalarField2D              b(num_nodes_x, num_nodes_y);

    ScalarEquation< CDS<double, 2> > single_phase(p, A, b, mesh.getDeltas());
    single_phase.setDeltaTime(dt); 
    single_phase.setGamma(Gamma);
    single_phase.setDirichlet(LEFT_WALL); 
    single_phase.setDirichlet(RIGHT_WALL);
    single_phase.setNeumann(TOP_WALL);
    single_phase.setNeumann(BOTTOM_WALL);
    single_phase.print();

    ofstream error_fp("./DataSingle2D/error"), residual_fp("./DataSingle2D/residual");
    int iteration = 0, Nerror = 0;
    double t = dt;
    double error = 1.0;
    double residual = 1.0;

    while (t <= Tmax) {
      std::cout << "\n ---- Iteration " << ++iteration 
		<< "\t time = " << t
		<< "\t Res = " << residual
		<< "\t Err = " << error;
      single_phase.calcCoefficients();
      Solver::TDMA2DX(single_phase, tolerance, tdma_iter, 1.0);
      residual = single_phase.calcResidual();
      //      error = single_phase.calcErrorL1();

      //      residual_fp << iteration << "\t" << residual << std::endl;	
      residual_fp << dt * iteration << "\t" << residual << std::endl;	
 
      single_phase.update();

      /* */
      error = 0.0;
      Nerror = 0;
      NumUtils::interpolateToNodes(pn, p);
      for(int i = pn.lbound(firstDim)+1; i <= pn.ubound(firstDim)-1; ++i)
	for(int j = pn.lbound(secondDim)+1; j <= pn.ubound(secondDim)-1; ++j) {

	  // L_2
	  error += (pn(i,j) - pa(i,j)) * (pn(i,j) - pa(i,j));
	  Nerror++;
	  // L_1
	  //	  if ( fabs( pn(i,j) - pa(i,j) ) > error )
	  //	    error = fabs (  pn(i,j) - pa(i,j) );
	}
      error = sqrt(error) / Nerror;
      //      error_fp << iteration << "\t" << error << std::endl;
      error_fp << dt * iteration << "\t" << error << std::endl;
      /* */

      t += dt;

      /*  */
      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(pn, p);
	InOut::writeToFile_DX(pn, iteration, "./DataSingle2D/pres.",dx, dy);	
      }
      /* */
    }

    std::cout << "\n\n NOTE: The data were written in the [./DataSingle2D] directory \n\n";
    return 0;
}




